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Abstract 

Spread options are a fundamental class of derivative contract written on multiple assets, 
and are widely used in a range of financial markets. There is a long history of approxima- 
tion methods for computing such products, but as yet there is no preferred approach that 
is accurate, efficient and flexible enough to apply in general models. The present paper 
introduces a new formula for general spread option pricing based on Fourier analysis of the 
spread option payoff function. Our detailed investigation proves the effectiveness of a fast 
Fourier transform implementation of this formula for the computation of prices. It is found 
to be easy to implement, stable, efficient and applicable in a wide variety of asset pricing 
models. 

Keywords: Spread options, multivariate spread options, jump-diffusions, fast Fourier trans- 
form, gamma function. 
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1 Introduction 

When Sjt,j = 1,2 are two asset price procsses, the basic spread option with maturity T and 
strike K > is the contract that pays (Sit — S2T — K)^ at time T. From the risk-neutral 
expectation formula, the time price of this option, assuming a constant interest rate r, must be 

Spr(5o; T, K) = e-''^ Es,[{S^t - S^t - K)+] . (1) 

The literature on applications of spread options is extensive, and is reviewed by Carmona 
and Durrleman yj who explore further applications of spread options beyond the case of equi- 
ties modelled by geometric Brownian motion, notably energy trading. For example, the "crack 
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spread" is the difference between the prices of crude oil and natural gas. "Spark spreads" refer 
to differences between the price of electricity and the price of fuel: options on spark spreads 
are widely used by power plant operators to optimize their revenue streams. Energy pricing re- 
quires models with mean reversion and jumps very different from geometric Brownian motion, 
and pricing spread options in such situations can be challenging. 

Closed formulas for ([T]) are known only for a limited set of cases. In the Bachelier stock 
model, St = {Su, 821) is an arithmetic Brownian motion, and in this case ([T]) has a Black- 
Scholes type formula for any T,K. In the special case K = when St is geometric Brownian 
motion, ([T]) is given by the Margrabe formula [|8l. 

In the important case where St is geometric Brownian motion and > 0, no explicit pricing 
formula is known. Therefore there is a long history of approximation methods. Numerical 
integration methods, typically Monte Carlo based, are often employed. Dempster and Hong [4ll 
introduced a numerical integration method based on double fast Fourier transforms (FFT) that 
is efficient for geometric Brownian motion and certain more general asset price processes. 

Another stream of research develops analytical methods applicable to log normal models 
that involve linear approximations of the nonlinear exercise boundary. Such methods are often 
very fast, but their accuracy is usually not easy to determine. Kirk [5J presented an analytical 
approximation that performs well in practice. Carmona and Durrleman [2J used a family of 
lower bounds for the spread option price to obtain an analytical approximation formula that 
apparently leads to acceptable accuracy. Moreover these results extend to approximate values 
for the Greeks. Deng, Li and Zhou [6] have developed a systematic analytical approximation 
scheme which is perhaps even more reliable. 

Next to having analytical option pricing formulas, the fastest option pricing engines by 
numerical integration are usually those based on the fast Fourier transform methods introduced 
by Carr and Madan [|3l. Their interest was in option pricing for geometric Levy process models 
like the variance gamma (VG) model, but their basic framework has been adapted to a host of 
other models where the characteristic function of the returns is known. 

The main purpose of the present paper is to give a new numerical integration method for 
computing spread options in two or higher dimensions via FFTs, applicable to any model for 
which the characteristic function of the joint return process is given analytically. Our method 
is completely different from prior approaches, and consequently has distinct advantages, no- 
tably the flexibility to include a wide range of desirable asset return models, with a competitive 
computational expense. 

The results we describe all stem from the following fundamental new result that gives a 
Fourier representation of the basic spread option payoff function P{xi, X2) = (e^'^ — e^'^ — 1)^- 
Note that this spread option payoff is without loss of generality: We can reduce the general case 

7^ to = 1 by using scaling and interchange of 5*1 and ^2. 

Theorem 1. For any real numbers e = (ei, £2) with £2 > and ei + 62 < —1 

PW = (2.)-// e'"'P(u)d\ Pin) = -jn-'u.) 

Here T{z) is the complex gamma function defined for '^e{z) > Q by the integral T{z) = 
e~^t^~^dt. We write u = {ui,U2) and x = (xi,X2) and the notation x' denotes the (un- 
conjugated) matrix transpose. 
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Using this result, whose proof is given in the appendix, we will find we can follow the 
logic of Carr and Madan to derive numerical algorithms for efficient computation of a variety 
of spread options and their Greeks. The basic strategy to compute ([T]) is to combine ([2]) with 
an explicit formula for the characteristic function of the bivariate random variable Xt = log St- 
For the remainder of this paper, we make a simplifying assumption that the dynamics of X is 
homogeneous, which is equivalent to the following factorization of the characteristic function 

ExoK^'^] = e™^o$(M; T), T) := Eo[e^"^^] . (3) 

Then, if St = e"^* is the price process for two traded assets, the spread option formula can be 
written 

Spr(5o;T) = e-'-^^K^iT - ^2T - 1)1 

= {27r)-^e-'^ f[ e'''''^^u;T)P{u)d\ . (4) 

J JR^+ie 

The Greeks are handled in exactly the same way. For example, the Delta := dSpr/dSio is 
obtained as a function of 5*0 by replacing $ in by d^/dSio- 

Explicit double Fourier integrals like this can be approximated numerically by a two dimen- 
sional FFT Such approximations involve both a truncation and discretization of the integral, 
and the two properties that determine their accuracy are the decay of the integrand of (|4]) in 
M-space, and the decay of the function Spr in x-space. The remaining issue of computing the 
gamma function is not a real difficulty. Fast and accurate computation of the complex gamma 
function in for example, Matlab, is based on the Lanczos approximation popularized by [9Q 

In this paper, we demonstrate how our method performs in three different models for spread 
options on two stocks, namely the geometric Brownian motion (GBM), a three factor stochastic 
volatility (SV) model and the variance gamma (VG) model. Section 2 provides the essential 
definitions of the three types of models, including explicit formulas for their bivariate charac- 
teristic functions. Section 3 discusses how the two dimensional FFT can be implemented for 
our problem, and gives a heuristic picture of how the accuracy and speed will depend on the 
choices made in the implementation. Section 4 gives the detailed results of the performance 
of the method in the three models. In that section, the accuracy of each model is compared to 
benchmark values computed by an independent method for a reference set of option prices. We 
also demonstrate that the computation of the spread option Greeks in such models is equally 
feasible. Section 5 extends all the above results to spread options on three or more assets. While 
in such cases the formulation is simple, the resulting higher dimensional FFTs are in practice 
much slower to compute. 

'According to these authors, computing the gamma function becomes "not much more difficult than other 
built-in functions that we take for granted, such as sin x or e^". 
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2 Three kinds of stock models 



2.1 The case of geometric Brownian motion 

In the two-asset Black-Scholes model, the vector St = (Su, 821) has components 

S,t = Sjoexp[{r-a^/2)t + a,W^],j = 1,2 

where cri,(T2 > and W^,W'^ are risk- neutral Brownian motions with constant correlation 
p, \p\ < 1. The joint characteristic function of Xt = (log Sit, log S'2t) as a function of m = 
{ui, U2) is of the form e*"^o$(M; T) with 

T) = exp[iu{rTe - o'^T 12)' - uT.u'T/2] (5) 

where e = (1, 1), S = [cr^, oxo^P', o"icr2p, cr|] and cr^ = diagS. Here we use implied matrix mul- 
tiplication, and the v! denotes the (unconjugated) matrix transpose. Substituting this expression 
into (|4]) yields the spread option formula 

j2 



Spr(S'o; T) = (27r)-2e-"^ / / e*"^o exp[w(rTe - a'^T 12)' - uLu'T /2]P{u)({' 



u . 



(6) 



As we discuss in Section 3, we recommend that this be computed numerically using the FFT. 



2.2 Three Factor Stochastic Volatihty Model 

The spread option problem in a three factor stochastic volatility model was given as an example 
by Dempster and Hong flU. Their model is defined by the SDE for X = (log Si, log 82)'- 

dXi = [{r-6i-cxl/2)dt + aiy^dW^] 
dX2 = [{r-62-(Tl/2)dt + a2VvdW^] 
dv = k{p — v)dt + av\/vdW^ 

where 

E[dW^dW^] = pdt 
EQldW^dW] = pidt 
EQldW^dW] = p2dt. 
As discussed in that paper, it has the joint characteristic function e'^'^^o^{u; T) where 
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-- [{(ylu\ + 0-2^^2 + 2p(TiCr2^ilM2) + i + Cr2'"2)] 

K - i{piaiUi + p2cr2U2)cru 
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2.3 Exponential Levy Models 



Many stock price models are of the form St = e'^* where Xt is a Levy process for which 
the characteristic function is explicitly known. We illustrate with example of the VG process 
introduced by [7], the three parameter process Yt with Levy characteristic triple (0, 0, u) where 
the Levy measure is z/(x) = A[e~"+^la;>o + e"-^'la;<o]/|x| for positive constants X,a±. The 
characteristic function of Yt is 



l + i 



1 



1 



u + 



u 



-At 



(7) 



To demonstrate the effects of correlation, we take a bivariate VG model driven by three in- 
dependent VG processes Yi, Y2, Y with common parameters a± and = = (1 — a)\, A^ = 
a\. The bivariate log return process Xt = log St is a mixture: 



X 



It 



XiQ + Yit + Yt] Xot — Xon + Y2t + Yt 



(8) 



Here a E [0, 1] leads to dependence between the two return processes, but leaves their marginal 
laws unchanged. An easy calculation leads to the bivariate characteristic function e*"^o$(M; T) 
with 



X 



1 + z 



l + t 










1 




a_ 


a+J 



1 



Ui 



ul 



Ul +U2) + ■ 



a_a+ 



,Ml + U2) 



a_aj 



l + i 



2 "I —a\t 



(9) 



a+J 



U2 + 



a_a_|_ 



-{l-a)\t 



3 Numerical Integration by Fast Fourier Transform 

To compute (|4]) in these models we approximate the double integral by a double sum over the 
lattice 

r = {u{k) = {U{ki),u{k2))\k = {ki,k2) G {0,...,Ar- 1}2}, u{k) = -U + k7] 

for appropriate choices of N,r],u := Nr]/2. For the FFT it is convenient to take to be a power 
of 2 and lattice spacing 77 such that truncation of the w-integrals to [—u, u] and discretization 
leads to an acceptable error. Finally, we choose initial values Xq = log 5*0 to lie on the reciprocal 
lattice with spacing 77* = 2Tc/Nr] = n/u: 

r = {xi£) = x(£2)) 1^ = (ii, £2) e {0, . . . , iV - 1}2}, xii) = -x + iri*, x = Nri*/2 . 

For any Sq = e"^" with Xq = x{£) E T* we then have the approximation 

2 -rT 

Spr(5o; T) ~ e^^"('=)+^^)^(')'$(«(fc) + te; T)P{u{k) + le) . (10) 

^ ' ki,k2=0 
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Now, as usual for the discrete FFT, as long as is even, 

iu{k)x{iy = in{ki + k2 + h + £2) + 2niki'/N (mod 27ii) . 
This leads to the double inverse discrete Fourier transform 

Spr(5o;T) ~ (-i)^i+^^e-^^ [^V e'^^^^)' ^ e'^^'^'' H {k) 

V / I., h„—n 
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fci,A;2=0 

-ly^+^^e-'^^ ^l^y e-^^(^)'[ifft2(i7)](^) (11) 
where 

H{k) = <!){u{k) + it- T)P{u{k) + it) . 

The selection of suitable values for e, and 77 is a somewhat subtle issue when implement- 
ing the above FFT approximation of (|4]). We now briefly discuss separately the truncation error 
and discretization error. The pure truncation error, defined by taking r] ^ 0, N ^ 00 keep- 
ing u = Nr]/2 fixed, can be made small if the integrand of (|4]) is small outside the square 
[—u + iti,u + iti] X [—u + it2,u + it2]. Corollary [4| proved in the Appendix, gives an 0(|u|~^) 
upper bound on P, while $(n) can generally been seen directly to have some n-decay. Typ- 
ically, with some caveats, if one picks u large enough so that |$| < 5i ^ 1 and has decay 
outside the square, then the truncation error will be less than 0{5i). 

The pure discretization error, defined by taking u ^ 00, N ^ 00 while keeping t] fixed, can 
be made small if e''"^"Spr, taken as a function of Xq G M^, has rapid decay in Xq. This is related 
to the smoothness of $(m) and the choice of t. The first two models are not very sensitive to e, 
but in the VG model the following conditions are needed to ensure that singularities in u space 
are avoided: 

— a+ < ei, 62, ei + 62 < a- . 
The error in the approximation ( [TT] ) with N = 00 and i] finite, is given by the formula 

Spr('')(Xo) - Spr(Xo) = e-'"^^'/''Spr(Xo + 2ix) . 

teZ2\{(0,0)} 

Typically, with some caveats, the norm of this can be made smaller than 0(52) ^ 1 for Xq E 
[—cx,cx]'^ if < c ^ 1 and e''"^oSpr(Xo) is 0(^2) and rapidly decaying outside the square 
[—X, x^. 

In summary, one expects that the combined truncation and discretization error will be 0{5i + 
^2) if M and r] = tt/x are each chosen as above. 

The FFT method can also be applied to the Greeks, enabling us to tackle hedging and other 
interesting problems. It is particularly efficient for the GBM model, where differentiation under 
the integral sign is always permissible. For instance, the FFT formula for vega (the sensitivity 
to a) takes the form: 

dSpr^ = (-l)^^+^^e-^^(|^ye-(^)'[ifft2(||)](£); 



dH{k) 
dai 



Mk) + It) [i^^^ + ^Hk) + ^t)'^ T/2 
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where = [2cri,0] and ^ = [2ai, pa2; p(T2,0]. Other Greeks including those of higher 
orders can be computed in similar fashion. This method needs to be used with care for the SV 
and VG models, since it is possible that differentiation leads to an integrand that decays slowly. 



Our numerical experiments were coded and implemented in Matlab version 7.6.0 on an Intel 
2.80 GHz machine running under Linux with 1 GB physical memory. If they were coded in C++ 
with similar algorithms, we should expect to see much faster performance. All the following 
experiments were run with ei = —3, 62 = 1. 

The strength of the FFT method is demonstrated by comparison with benchmarks. Based 
on a selection of {SIq, S^q, K^},i G 1, 2, 3 ... n, the objective function we measure is defined 
as 



where M* and 5* are FFT computed prices and benchmark prices for the ith combination. 
We chose log 5*10 = ^, log 5*20 = -f + j G 1,2,3,4,5,6 and fixed the strike K = 1, 
leading to 36 prices contributing to the objective function. These choices cover a wide range 
of moneyness, from deep out-of-the-money to deep in-the-money. Since these combinations all 
lie on lattices T* corresponding to = 2" and m/10 = 2™ for integers n, m, all 36 prices M* 
can be computed simultaneously with one FFT. 

Figure [T] shows how the FFT method performs in the 2-dimensional Geometric Brownian 
motion model for different choices of N and u. Since the two factors are bivariate normal, 
benchmark prices can be calculated to high accuracy by one dimensional integrations. In Figure 
[T]we can clearly see the effects of both truncation errors and discretization errors. For a fixed 
u, the objective function decreases when increases. The n = 20 curve flattens out near 10~^ 
due to its truncation error of that magnitude. In turn, we can quantify its discretization errors 
with respect to A^ by subtracting the truncation error from the total error. The flattening of the 
other three curves near 10~^^ should be attributed to Matlab round-off errors: because of the 
rapid decrease of the characteristic function $, their truncation error is negligible. For a fixed 
A^, increasing u brings two effects: reducing truncation error and enlarging discretization error. 
These effects are well demonstrated in figure [T| 

For the stochastic volatility model, no analytical or numerical method we know is consis- 
tently more accurate to serve as benchmark. Instead, we used Monte Carlo simulation to achieve 
a moderate degree of accuracy. Each benchmark price was created by 1, 000, 000 simulations, 
each consisting of 2000 time steps. No variance reduction was employed. The resulting ob- 
jective function is shown in Figure [Ij The objective function only decreases to near 4 x 10^^, 
reflecting the inaccuracy of the benchmark. A comparable graph (not shown), using benchmark 
prices computed with the FFT method with A^ = 2^^ and u = 80, showed similar behaviour to 
Figure [T] and is consistent with accuracies at the level of round-off. The computational cost to 
reduce the Monte Carlo simulation error to a comparable level would be prohibitive. 

The VG model also lacks a reliable benchmark, so again we used Monte Carlo simulation. 
Since the VG process can be simulated as a Brownian motion time-changed by a gamma pro- 



4 Numerical Results 
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cess, MC is here much more efficient than for the SV model, and we were able to run 10^ 
simulations with antithetic variates for each price. The objective function is shown in Figure 
[3] The truncation error for m = 20 is about 2 x 10~^. The other three curves flatten out 
near 7 x 10~^, which is presumably the accuracy of the benchmark. A comparable graph (not 
shown), using benchmark prices computed with the FFT method with = 2^^ and u = 80, 
showed similar behaviour to Figure [1} and is consistent with the FFT method being capable of 
producing accuracies at the level of round-off. 

The strength of the FFT method is further illustrated by the computation of individual prices 
shown in Tables [T] to [3} The columns labeled "Analytic" and "Simulation" refer to the type of 
benchmark used. One can observe that an FFT with N = 256 is capable of producing very high 
accuracy in all three models. In the GBM model, the errors are simply due to round-off. The 
discrepancies in the VG model are small too, less than 1 basis point. The apparent discrepancies 
in the SV model we attribute to the inaccurate benchmark, rather than the FFT method itself. 

The FFT computes in a single iteration an x panel of prices spread corresponding 
to initial values Sio = e^'i»+^i^*, ^20 = e''^°+^^'^\ K = 1, (£1,4) G {0,...,A^ - 1}\ If 
the desired selection of {Sio, S20, K} fits into this panel of prices, or its scaling, a single FFT 
suffices. If not, then one has to match (xio, 0:20) with each combination, and run several FFTs, 
with a consequent increase in computation time. However, we have found that an interpolation 
technique is very accurate for practical purposes. For instance, prices for multiple strikes with 
the same and ^20 are approximated by a polynomial fit along the diagonal of the price panel: 
Spr(5'o; Ki) = Krspread{l, 1), Spr(5o; Kle-''^*) = Kie-"^' ■spread{2, 2), Spr(5o; J^ie-^'?*) = 
Kie^'^'^* ■ spread{3, 3) . . . . The results of this technique are recorded in Table|2]and Tablejsjin 
the column "Interpolation". We can see this technique generates very competitive results and 
moreover, saves computational resource. 

Finally, we computed first order Greeks using the method described at the beginning of 
Section|3]and compared them with finite differences. As seen in Table|4} the two methods come 
up with very consistent results. The Greeks of our at-the-money spread option exhibit some 
resemblance to the at-the-money European put/call option. The delta of 5*1 is close to the delta 
of the call option, which is about 0.5. And the delta of S2 is close to the delta of the put option, 
which is also about 0.5. The time premium of the spread option is positive. The option price 
is much more sensitive to 5*1 volatility than to 5*2 volatility. It is an important feature that the 
option price is negatively correlated with the underlying correlation: Intuitively speaking, if the 
two underlyings are strongly correlated, their co-movements diminish the probability that Sit 
develops a wide spread over S2t- This result is consistent with observations made by 0. 

Since the FFT method naturally generates a panel of prices, and interpolation can be imple- 
mented accurately with negligible additional computational cost, it is appropriate to measure 
the efficiency of the method by timing the computation of a panel of prices. Such computing 
times are shown in Table [5| For the FFT method, the main computational cost comes from 
the calculation of the matrix H in ( fTT] ) and the subsequent FFT of H. We see that the GBM 
model is noticeably faster than the SV and VG models: This is due to a recursive method used 
to calculate the H matrix entries of the GBM model, which is not applicable for the SV and VG 
models. The number of calculations for H is of order A^^ which for large A^ exceeds the A^ log A^ 
of the FFT of H, and thus the advantage of this efficient algorithm for GBM is magnified as A^ 
increases. However, our FFT method is still very fast for the SV and VG models and is able to 



8 



generate a large panel of prices within a couple of seconds. 



5 High Dimensional Spread Options 

The ideas of section 2 turn out to extend naturally to basket options with payoffs {St — Sit — 
■ ■ ■ — Smt — l)"*" for M > 2. The analysis is based on a simple result proved in the Appendix: 

Proposition 2. Let z E M. and u = {ui, . . . , um)' £ with '^m{um) > Ofor all m < M. 
Then 



M 



m=l ^( ^ Sm=l "^"i) 

As before we write Xm = log S'm, x = log S and seek to compute for u e C^'^, u E C 

AI 



(13) 



P(u,u) 



JRM+l 



We introduce the factor 1 = 8{e^ — X]m=i G^"^)G^dz and interchange the z integral with the 
X integrals. Then using Proposition [2] one finds 



Piu.u) 



„ M 



-tux(^x _ _ 1 



„ M 



dxdz 



dxdz 



1^( ^ Sm=l 

When we can apply Proposition 1, we obtain the result: 

Proposition 3. For any real numbers e, e = (ei, . . . , eu) with em > for all m < M and 



M 



m=l 



i2n 



where 



.-(M+l) 



r(m + 1) 



P{u,u)d^'udu (14) 



MA-r+i+i(e,g) 



(15) 
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6 Conclusion 



This paper presents a fundamentally new approach to the valuation of spread options, an impor- 
tant class of financial contracts. The method is based on a newly discovered explicit formula 
for the Fourier transform of the spread option payoff in terms of the gamma function. 

This mathematical result leads to simple and transparent algorithms for computing spread 
options in all dimensions. The powerful tool of the Fast Fourier Transform provides an accurate 
and efficient implementation of the fundamental result. The difficulties and pitfalls of the FFT, 
of which there are admittedly several, are by now well understood, and thus the reliability and 
stability properties of our method are rather clear. 

Many important processes in finance, particularly affine models and Levy jump models, 
have well known explicit characteristic functions, and can be included in the method with little 
difficulty. Thus the method can be easily applied to important problems arising in energy and 
commodity markets. 

Finally, the Greeks can be systematically evaluated for such models, with similar perfor- 
mance and little extra work. 

A Proof of Theorem |T] and Proposition |2| 

Proof of Theorem [l| Suppose £2 > 0, ci + £2 < —1. One can then verify either directly or 
from the argument that follows that e^'^P{x), e = (ei, 62) is in L^(]R^). Therefore, application 
of the Fourier inversion theorem to e'''^P{x),e = (ei, 62) implies that 



P{x) = (27r)" 



(16) 



where 

g(u) = [[ e-*"-^P(x)rf2a: . 
J Jr^ 

By restricting to the domain {x : xi > 0, e^^ < e"^^ — 1} we have 

»log{e^l-l) 



9[u) 
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-iu2 1 — iU2 
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The change of variables 2; = e then leads to 
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The beta function 
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is defined for any complex a, b with 3fte(a), 3?e(6) > by 

B{a,b) = [ z^-^l - zf-^dz 







From this, and the property T{z) = {z — l)T{z — 1) follow the formulas 

^ r(z(m + n2)-l)r(-m2 + 2) ^ r(z(Mi + n2) - i)r(-m2) 
^^""^ (l-m2)(-m2)r(mi + l) r(mi + l) ' ^ ^ 

□ 

The above derivation also leads to the following bound on P. 
Corollary 4. Fix t2 = e, ei = — 1 — 2tfor some e > 0. Then 

r(6)r(2 + 6) ^ 

r(2 + 2e) Q(|m|V5)i/2 
w/zere = + e^){z + (1 + e)^). 

Proof: First note that for zi,Z2 E C, \B{zi,Z2)\ < B{^e{zi),^e{z2)). Then ([T7]) and a 
symmetric formula with U2 —I — ui — U2 leads to the upper bound 

1 



|P(mi -i(e + 1),U2 + ie)\ < B{e,2 + e) min 



Q{\U2\)'Q{\UI+U2\) 



But since Q is monotonic and |m| < max (|m2|, \ui + W2I) for all m G M^, the required result 
follows. □ 

Proof of Proposition |2| To simplify the proof, we make the change of variables p = and 
Qm = e"^™ . It becomes to prove that 

/ ps{p - E n c"-'^''^ = l:"'=\tM '^"^> -(^-- . (19) 

We prove the proposition by induction. The above equation trivially holds when M = 1. 
Now suppose it holds for M = N. Then for M = iV + 1 

„ N N 

LHS = / - g^+i - 5^ g Jg^+T^"' n """'^"^^'^ 

\ m=l 



11^=1 r(-^Mm) r P 



r(-^Lm=i^™) >/0 P-QN+l 

The above step makes use of the result when M = N. The proof is then complete when 
one notices that the ^at+i integral is simply p"Hz.m=i ""i-i multiplied by a beta function with 
parameters —i{J2m=i'^m) and —iuN+i- 

□ 
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Error Plot for Spread Option (GBM) 




Figure 1 : This graph shows the objective function Err for the numerical computation of the 
GBM spread option versus the benchmark. Errors are plotted against the grid size for different 
choices of u. The parameter values are taken from [4J: r = 0.1, T = 1.0, p = 0.5, (5i = 
0.05, (Ti = 0.2, 52 = 0.05, era = 0.1. 



Table 1 : Prices for the two-factor GBM model of H for different choices of A^. The parameter 
values are the same as Figure [1] except we fix Sio = 100, S20 = 96,u = 40. 



Strike K 


Analytic 


64 


128 


256 


512 


0.4 


8.312461 


8.206666 


8.312331 


8.312461 


8.312461 


0.8 


8.114994 


8.009643 


8.114864 


8.114994 


8.114994 


1.2 


7.920820 


7.815913 


7.920691 


7.920820 


7.920820 


1.6 


7.729932 


7.625469 


7.729804 


7.729932 


7.729932 


2.0 


7.542324 


7.438304 


7.542196 


7.542324 


7.542324 


2.4 


7.357984 


7.254408 


7.357857 


7.357984 


7.357984 


2.8 


7.176902 


7.073770 


7.176775 


7.176902 


7.176902 


3.2 


6.999065 


6.896377 


6.998939 


6.999065 


6.999065 


3.6 


6.824458 


6.722213 


6.824332 


6.824458 


6.824458 


4.0 


6.653065 


6.551264 


6.652940 


6.653065 


6.653065 
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Error Plot for Spread Option (SV) 
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Figure 2: This graph shows the objective function Err for the numerical computation of the SV 
spread option versus the benchmark Monte Carlo simulation values computed with 1, 000, 000 
simulations each of which consists of 2000 time steps. The parameter values are taken from 
S: r = 0.1, T = 1.0, p = 0.5, 5i = 0.05, ai = 1.0, pi = -0.5,^2 = 0.05,^2 = 0.5, pa = 
0.25, vo = 0.04, K = 1.0, p = 0.04, = 0.05. 
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Error Plot for Spread Option (VG) 



10' 




Figure 3: This graph shows the objective function Err for the numerical computation of the 
VG spread option versus the benchmark values computed with Monte Carlo simulation, where 
10^ simulations are performed for each price. Errors are plotted against the grid size for five 
different choices of u. The parameters are: r = 0.1, T = 1.0, p = 0.5, a+ = 20.4499, a_ = 
24.4499, a = 0.4, A = 10 
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Table 2: Prices for the 3 factor SV model of H for different choices of N. The parameter 
values are the same as Figure [2] except we fix 5*10 = 100, 5*20 = 96, m = 40. The interpolation 
is based on a matrix with discretization of = 256 and a polynomial with degree of 8. 



oLTlKc JV 




1 9R 

iZO 






illLcrpOldLlOIl 


OllIlUldLlOIl 


2.0 


0.996467 


n Z A A OCT 

7.544853 


7.548502 


7.548502 


n z AO cr\'^ 

7.548502 


7.545784 


2.2 


6.902676 


7.449895 


7.453536 


7.453536 


7.453536 


7.435475 


2.4 


6.809696 


7.355748 


7.359381 


7.359381 


7.359381 


7.360184 


2.6 


6.717527 


7.262411 


7.266036 


7.266037 


7.266036 


7.276051 


2.8 


6.626167 


7.169883 


7.173501 


7.173501 


7.173501 


7.193546 


3.0 


6.535616 


7.078165 


7.081775 


7.081775 


7.081775 


7.077421 


3.2 


6.445873 


6.987254 


6.990856 


6.990857 


6.990856 


6.983995 


3.4 


6.356936 


6.897150 


6.900745 


6.900745 


6.900745 


6.913994 


3.6 


6.268806 


6.807853 


6.811439 


6.811440 


6.811439 


6.818526 


3.8 


6.181481 


6.719360 


6.722939 


6.722939 


6.722939 


6.735112 


4.0 


6.094959 


6.631670 


6.635241 


6.635242 


6.635241 


6.621617 



Table 3: Prices for the VG model of for different choices of A^. The parameter values are the 
same as Figure |3] except we fix Siq = 100, 5*20 = 96, u = 40. The interpolation is based on a 
matrix with discretization of = 256 and a polynomial with degree of 8. For the last column, 
10^ simulations are used for each price. 



Strike K 


64 


128 


256 


512 


Interpolation 


Simulation 


2.0 


9.157674 


9.723691 


9.727458 


9.727458 


9.727458 


9.727546 


2.2 


9.061487 


9.626247 


9.630006 


9.630006 


9.630006 


9.629444 


2.4 


8.965876 


9.529448 


9.533200 


9.533200 


9.533200 


9.533137 


2.6 


8.870896 


9.433296 


9.437040 


9.437040 


9.437040 


9.437323 


2.8 


8.776560 


9.337792 


9.341527 


9.341528 


9.341527 


9.341547 


3.0 


8.682870 


9.242934 


9.246662 


9.246662 


9.246662 


9.246273 


3.2 


8.589828 


9.148725 


9.152445 


9.152445 


9.152445 


9.152956 


3.4 


8.497433 


9.055163 


9.058875 


9.058875 


9.058875 


9.058997 


3.6 


8.405687 


8.962250 


8.965954 


8.965954 


8.965954 


8.966410 


3.8 


8.314590 


8.869984 


8.873681 


8.873681 


8.873681 


8.874129 


4.0 


8.224140 


8.778368 


8.782057 


8.782057 


8.782057 


8.781863 



Table 4: The Greeks for the GBM model compared between the FFT method and the finite 
difference method. The FFT method uses N = 2^° and u = 40. The finite difference uses a 
two-point central formula, in which the displacement is ±1%. Other parameters are the same 
as Table [T] except that we fix the strike K = A.O to make the option at-the-money. 





Delta(Sl) 


Delta(S2) 


Theta 


Vega((Ti) 


Vega((T2) 


dSpr/dp 


FD 


0.512648 


-0.447127 


3.023823 


33.114315 


-0.798959 


-4.193749 


FFT 


0.512705 


-0.447079 


3.023777 


33.114834 


-0.798972 


-4.193728 
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Table 5: Computing time of FFT for a panel of prices. 



Discretization 


GBM 


SV 


VG 


64 


0.091647 


0.083326 


0.109537 


128 


0.099994 


0.120412 


0.139276 


256 


0.126687 


0.234024 


0.220364 


512 


0.240938 


0.711395 


0.621074 


1024 


0.609860 


2.628901 


2.208770 


2048 


2.261325 


10.243228 


8.695122 
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